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SUMMARY 


A  numerical  procedure  for  solving  the  problem  of 
steady  supersonic  inviscid  flow  around  smooth  conical  bodies  is 
presented.  Results  are  obtained  by  solving  the  elliptic  partial 
differential  equations  that  define  the  conical  flow  between  the 
body  and  the  shock.  Results  are  given  for  circular  cones  up  to 
moderately  high  relative  incidences,  including  some  cases  for 
incidences  beyond  a  critical  value  at  which  the  entropy  sin¬ 
gularity  moves  from  the  surface. 

Also  presented  are  a  few  results  for  elliptic  cones  at 
zero  and  non-zero  incidence,  as  w'ell  as  results  for  another 
conical  body  whose  cross  section  is  defined  by  a  fourth  order 
even  cosine  Fourier  series. 

The  applicability  of  the  method  can  be  limited  by  the 
entropy  singularity  moving  too  far  away  from  the  surface,  by 
the  flow  field  containing  regions  of  locally  conically  supersonic 
flow,  or  by  the  shock  wave  approaching  very  close  to  the  Mach 
wave. 


Comparison  of  results  shows  excellent  agreement  with 
other  theoretical  methods  and  also  with  experimental  results. 
The  method  is  efficient  in  computer  time. 
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1.0  INTRODUCTION 

The  first  attempts  to  obtain  solutions  for  the  steady  inviscid  supersonic  flow  field  about 
conical  bodies  were  made  by  several  authors  in  the  early  1930’s.  These  authors  considered  a  circular 
cone  at  zero  incidence  to  the  free  stream.  The  equations  defining  the  flow  were  reduced  to  ordinary 
non-linear  differential  equations  to  be  satisfied  between  the  shock  and  the  body,  thus  giving  a  two- 
point  boundary  value  problem.  These  equations  were  numerically  integrated  by  Kopal  in  1947 
(Ref.  1)  using  the  equations  as  derived  by  Taylor  and  Macoll  in  1932.  Later,  Sims  computed  these 
zero  order  solutions  again  and  presented  the  results  in  a  convenient  form  in  Reference  2. 

Using  a  small  perturbation  method,  Stone  in  1948  (Ref.  3)  produced  a  first  order  solution 
applicable  to  circular  cones  at  small  incidence.  Kopal  made  the  original  numerical  integration  of  the 
resulting  equations,  and  Sims  (Ref.  4)  later  repeated  and  extended  the  calculations.  A  second  order 
theory  for  a  circular  cone  was  also  formulated  by  Stone  (Ref.  5),  and  Kopal  in  1949  again  provided 
numerical  results.  These  second  order  tabulations  of  Kopal  (Ref.  6)  are  of  limited  range  and  in¬ 
convenient  to  use. 


These  first  and  second  order  solutions  contain  uncertainties,  as  pointed  out  by  Ferri  (Ref. 
7),  owing  to  the  singularity  in  entropy  not  being  included.  Ferri  observed  that  the  body  must  be  a 
surface  of  constant  entropy  and  that  all  streamlines  must  eventually  converge  at  the  leeward  gen¬ 
erator  of  the  cone,  thus  resulting  in  an  entropy  singularity  at  this  line.  Ferri  also  pointed  out  that 
there  is  a  layer  very  close  to  the  surface,  the  vortical. layer,  in  which  the  entropy  changes  rapidly  in 
a  direction  normal  to  the  surface. 

Melnik,  amongst  others,  studied  in  detail  the  vortical  layer  and  the  position  of  the  entropy 
singularities  (Ref.  8  and  9)  for  a  circular  cone  and  a  delta  whig  and  was  able  to  show  that,  hi  the 
case  of  the  circular  cone,  the  nodal  type  singularity  remained  at  the  leeward  generator  until  an  in¬ 
cidence  was  reached  when  the  transverse  pressure  gradient  became  adverse.  Beyond  this  incidence 
it  was  shown  that  the  nodal  singularity  might  follow  the  point  of  minimum  pressure  or  might  stay 
at  the  leeward  generator  until  a  certain  incidence  was  reached,  at  which  time  the  singularity  would 
lift  off  the  cone  surface.  It  is  shown  in  this  Report  that  the  latter  is  the  correct  behaviour. 


Several  numerical  techniques  have  been  developed  in  recent  years  to  solve  the  problem, 
particularly  for  cones  that  are  circular  or  elliptical  in  cross  section.  Stocker  and  Mauger  (Ref.  10), 
amongst  others,  solved  the  set  of  elliptic  partial  differential  equations  that  apply  between  the  shock 
wave  and  the  body,  assuming  a  shock  shape  and  integrating  in  towards  a  body  that  is  defined  by 
the  envelope  of  the  streamlines  obtained  in  the  integration  procedure.  The  shock  shape  was  then 
modified  empirically  to  obtain  a  given  body  shape.  Results  were  presented  for  the  circular  cone  with 
a  semi-apex  angle  of  20°  (at  incidences  of  5°  and  10°)  for  a  Mach  number  of  3.53.  The  results  for  the 
5°  incidence  agree  well  with  experiment,  while  only  fair  agreement  was  obtained  at  10°  incidence.  In 
fact,  difficulties  were  encountered  in  the  iteration  procedure  for  the  10°  incidence  case  and  it  was 
not  possible  to  obtain  a  solution  for  an  incidence  of  15°.  At  this  larger  incidence,  the  computed  body 
seemed  to  have  a  “bump”  located  at  and  near  the  leeward  generator,  which  at  that  time  was  thought 
to  be  accounted  for  by  the  singularity  leaving  the  body.  Stocker  and  Mauger’s  method  can' also,  in 
principle,  be  applied  to  other  bodies '  besides  circular  cones,  and  a  computation  was  made  for  an 
elliptic  cone  at  zero  incidence,  the  results  of  which  gave  reasonable  agreement  with  experiment. 
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A  different  numerical  method  of  solving  the  problem  for  any  conical  body  was  proposed  by 
Babenko  et  al.  (Ref.  11).  In  this  method  the  full  three-dimensional  flow  field  equations  are  considered 
and,  starting  from  an  estimated  shock  shape  with  the  given  body  set  at  the  required  incidence,  the 
equations  are  integrated  step  by  step  downstream  until  a  condition  of  conicity  is  reached.  When 
this  condition  is  reached,  the  full  flow  field  solution  is  then  available. 

Reference  11  gives  tabulated  values  of  the  flow  field  quantities  for  circular  cones  at  incidence. 
Some  results  are  given  for  Mach  numbers  of  2,  3,  4,  5,  6,  and  7,  cone  semi-angles  of  from  10°  to  45° 
in  5°  steps,  and  for  relative  incidences  up  to  0.8.  Most  of  the  Mach  number  4  and  6  results  were 
obtained  by  interpolation  from  the  results  at  the  other  Mach  numbers.  These  numerioal  results 
show  excellent  agreement  with  experiment. 

Gonidou  (Ref.  12)  has  used  the  method  of  Babenko  et  al.  to  obtain  solutions  for  a  circular 
cone  at  relative  incidence,  a/0o,  up  to  about  1.2  where  the  entropy  singularity  has  lifted  from  the 
surface  of  the  cone  at  the  leeward  generator.  Solutions  for  elliptic  cones  of  fairly  high  eccentricity 
(~3)  were  also  obtained  by  Gonidou.  Computer  times  and  some  other  program  details  for  the 
Babenko  method  are  given  in  his  paper.  About  500  downstream  steps  are  required  to  reach  a  condi¬ 
tion  of  conicity,  and  about  Y  hour  is  required  for  a  typical  solution  on  a  CDC  3600  computer  for  a 
mesh  size  of  10  x  16  in  the  radial  and  circumferential  directions  respectively. 

Moretti  (Ref.  13)  has  used  a  similar  approach  to  that  of  Babenko  et  al.  Again,  the  flow 
field  solution  is  obtained  by  marching  step '  by  step  downstream  until  a  conicity  condition  is  suffi¬ 
ciently  well  satisfied.  The  method  differs  from  Babenko’s  method  in  the  details  of  numerical  analysis 
and  in  his  use  of  a  characteristics  method  on  the  shock  and  on  the  body.  Moretti’s  method  requires 
about  400  downstream  steps  and  the  computer  time  is  typically  }/%  hour  on  an  IBM  360/50  computer. 
This  time  is  applicable  to  a  mesh  increment  that  is  about  twice  the  size  of  that  used  by  Gonidou  in 
the  time  quoted  above.  To  get  a  true  comparison  of  computer  times,  it  should  be  noted  that  a  CDC 
3600  computer  operates  at  about  twice  the  speed  of  an  IBM  360/50  computer. 

The  purpose  of  investigating  a  further  numerical  technique,  as  given  in  this  Report,  was 
not  only  to  try  to  find  a  more  efficient  technique  than  those  given  previously,  but  also  to  investigate 
a  method  that,  in  principle,  is  capable  of  solving  non-linear  elliptic  partial  differential  equations;  in 
itself  a  difficult  problem  in  numerical  analysis.  The  present  method  uses  the  condition  of  conicity 
to  reduce  the  problem  to  a  set  of  elliptic  non-linear  partial  differential  equations  in  two  independent 
variables.  A  transformation  of  co-ordinates  is  used,  as  in  the  methods  of  Babenko  and  Moretti,  to 
fix  the  boundaries  between  which  the  elliptic  equations  are  to  be  satisfied.  This  transformation  also 
has  the  effect  of  including  the  body  shape  in  the  coefficients  of  the  partial  differential  equations  and 
in  the  boundary  conditions,  so  that  the  same  method  can  be  used  for  general  conical  body  shapes . 
simply  by  changing  a  few  program  statements  to  redefine  the  equation  of  the  body.  In  fact,  the 
method  is,  in  most  cases,  only  limited  by  locally  supersonic  cross-flow  conditions,  or  by  the  entropy 
singularity  moving  too  far  away  from  the  surface,  or  by  the  shock  approaching  very  close  to  the 
Mach  wave. 

At  the  present  time  the  method  has  been  used  successfully  for  circular  cones  and  for  bodies 
that  can  be  obtained  by  successive  perturbations  of  a  circular  cone  and  that  do  not  have  curvatures 
that  are  too  large.  The  examples  given  here  are  for  circular  cones  at  incidence,  elliptic  cones,  and  a 
body  whose  cross-sectional  shape  is  represented  by  a  fourth  order  even  cosine  Fourier  series. 


t 
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The  method  is  efficient  in  computer  time  compared  with  other  fully  numerical  techniques 
and  one  solution  takes  from  about  Yi  minute'  to  3  minutes  on  an  IBM  360/50  computer  for  the 
circular  cone  at  incidence  —  the  time  increasing  as  .the  incidence  increases. 

2.0  EQUATIONS  OF  MOTION  AND  THE  BOUNDARY  CONDITIONS 


2.1  The  Equations  of  Motion  and  the  Boundary  Conditions  for  a  General  Body 

The  co-ordinate  system  and  equations  of  motion  are  Written  in  a  notation  similar  to  that 
used  in  Reference  11. 


Let  (z,  r,  9)  be  a  cylindrical  co-ordinate  system  as  shown  in  Figure  1.  Then  the .  equations 
of  continuity,  momentum,  and  energy  for  an  inviscid,  non-heat  conducting  gas  can  be  written  in 
matrix  form  in  this  co-ordinate  system  as  follows 


where 
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In  the  above  matrices  and  vectors  (u,  v,  w)  are  the  velocity  components  in  the  (z,  r,  9) 
directions  respectively,  p  is  the  pressure,  p  the  density,  and  c‘-  the  square  of  the  local  speed  of  sound. 


Boundary  conditions  are  given  at  the  body  and  at  the  shock.  The  boundary  condition  on 
the  body  is  that  the  normal  velocity  should  be  zero.  Thus,  if  the  equation  of  the  body  is  r  =  g  (z,  6) 
then  the  boundary  condition  on  this  body  is 


dg 
u  t5- 
dz 


V  + 


W  dg 


g  96 


=  0 


(2) 
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Velocity  components,  pressure,  and  density  behind  the  unknown  shock  can  be  found  in 
terms  of  the  shock  equation  by  applying  the  Rankine-Hugoniot  relations  across  the  discontinuity, 
giving  the  following  equations 

pV »  =  pm  V,m  ■  . 

P  +  Pco  Vn<X,  V»  =  Pco  +  Pco 

h  +  y2Vn  =h  m+  y2  v„l  (3) 


■  df  ,  df 

u  +  v^  =  u00  +Voo- 


v  af  df 

fdd  +  w=  T~dd  +  w- 


where  Y  is  the  velocity  component  nonnal  to  and  behind  the  shock,  and  is  given  by 


Vn  = 


V  + 


w  df 

r  so 


*y-  +  ii  si\ 

az )  +  \f  as) 


h  is  the  static  enthalpy  (=  S  for  a  calorically  perfect  gas)  and  the  equation  of  the  shock  is  given 
by  r  =  f  (z,  8).  The  subscript  “  refers  to  values  in  front  of  the  shock  wave,  hence 

Uco  =  Yc  COS  «»  Voo  =  -  Yo  sin  «  C0S  6>  Woo  =  Yo  sin  «  sin  0 


and 


where  Yo  is  the  free  stream  velocity. 


The  present  Report  deals  with  situations  that  are  conical,  hence  the  number  of  independent 
variables  can  be  reduced  from  three  to  two.  A  suitable  transformation  to  carry  out  this  reduction  is 
made  in  the  following  paragraph. 


2.2  Transformation  of  the  Independent  Variables 

The  equations  given  in  the  previous  paragraph  apply  to  any  general  body  shape  r  = 
g  (z,  8)  and  a  shock  wave  given  by  r  =  f  (z,  6).  However,  we  now  consider  a  situation  in  which  the 
body  and  flow  characteristics  are  conical  so  that  the  number  of  independent  variables  can  be  reduced 
from  three  to  two. 

Firstly  it  is  seen  that,  since  the  z-axis  is  along  the  axis  of  the  cone  with  the  origin  at  the 
apex,  the  equation  of  the  conical  body  can  be  written 


r  =  z  G(0 ) 
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where  G  (0)  is  a  function  of  d  only,  and  similarly  the  equation  of  the  conical  shock  wave  can  be  written 

r  =  z  F(0) 

where  F(0)  is  a  function  of'  0  only. 

Now  the  independent  variables  (z,  r,  0)  are  transformed  by  transformations  involving  the 
equations  of  the  body  and  the  shock  as  follows 

x  =  z 


,  _ _ r  -  z  G  (0) 

5  z[F  (0)  -  G  (0)  ] 

<t>  =  e 

This  transformation  has  the  effect  of  fixing  the  radial  as  well  as  the  circumferential  boundaries 
between  which  the  partial  differential  equations  are  to  be  satisfied,  since  now  the  body  r  =  z  G(0) 
corresponds  to  the  surface  £  =  0,  while  the  shock  r  =  z  F(0)  corresponds  to  the  surface  £  =  1. 


The  above  transformations  change  the  equations  of  motion  (1)  to  the  form 

Af  +  Bf  +  0f  +  D-° 

where  A  —  A' 

B  =  G  A'  +  £r  B'  +  £fl  C' 

C  =  C’ 

D  =  D' 

and  the  subscripts  denote  partial  differentiation  so  that  £z,  £r,  and  £  can  be  evaluated  as 


£  r 


z  (F  -  G) 

£z  =  -  £r  [G  +  (F  -  G)£]  , 

k  =  ~  z  ^  [Ge  +  -  Ge>  €1 


(4) 


In  the  new  co-ordinate  system  (x,  £,  4>)  the  region  in  which  the  equations  are  to  be  satisfied  is 

x  >  0 


0  <  £  <  1 
0  <  <p  <2  ir 

However,  since  the  flow  is  conical,  quantities  are  constant  for  all  x  and  for  a  fixed  £  and  <j>,  thus 
c)X 

=  0.  The  problem  is  then  reduced  to  two  dimensions  in  the  independent  variables  £  and  4> 
and  we  may  consider  the  equations  at  a  unit  distance  (x  =  z  —  1)  along  the  body. 


—  6  — 


Thus  the  problem  is  reduced  to  that  of  solving  the  elliptic  non-linear  partial  differential 


equations 


B  +  C  ^  +  D  =  0 

Ol ;  d(j) 


(5) 


subject  to  the  boundary  conditions  given  in  the  next  paragraph. 


2.3  The  Boundary  Conditions  for  a  Conical  Body 

In  the  transformed  co-ordinate  system  the  boundary  conditions  become,  from  equation  (2), 


at  £  =  0 


,  w  dG  A 
uG  -v+G'o*  =0 


(6) 


and  at  £  =  1,  from  equations  (3) 

pVIi  —  pm  Vi, 


P  +  Pro  Vna)  Vn  —  P^  -f-p^  Vnm 

h  +  hv;  =  h^  +  y2Xm. 


u+vF  =  ueo  +  vctiF 


v  t)F  ,  _  Voo  aF 

F  d<i>  +  w  F  d<t>  +  w“ 


(V) 


where 


„  ,  W  dF 

uF-v+f  -- 

y  _ _ E _ 9(j> 

and  the  equations  of  the  body  and  shock  have  transformed  to 

r  =  GM 

and  r  =  F(</>) 

respectively. 

Further  boundary  conditions  can  be  applied  if  the  body  and  flow  characteristics  are  sym¬ 
metrical  about  the  axis  6(  =  cj>)  =  0,  r  as  is  often  the  case.  In  this  case  boundary  conditions  at  </>  = 
0,  7r  are 

du  _  3v  _  dp  _  dp  _  —  0 

dcj)  dcj)  dcj)  dcj>  W 

and  the  problem  is  reduced  to  solving  the  equations  in  the  region 

0  <  ?  <  1 
0  IS  <j)  <  7T 


In  the  next  Section  a  summary  of  the  method  for  solving  the  set  (5),  subject  to  the  given 
boundary  conditions,  is  presented;  a  more  detailed  description  is  presented  in  Appendix  A, 


3.0  METHOD  OF  SOLUTION 


We  observe  that  if  the  function  F (<£)  defining  the  shock  shape  is  known,  the  problem  is 
solved  completely,  since  then  the  equations  (7)  can  be  solved  to  give  u,  v,  w,  p,  p  at  £  =  1  and  equa- 

tions  (5)  can  be  integrated  from  £  =  1  to  £  =  0  by  replacing  the  elements  of  .  by  difference  formulae 

()(j) 

of  the  type(A-l)  (see  Appendix  A).  This  makes  the  set  (5)  into  ordinary  differential  equations  and 
they  can  be  integrated  by  standard  techniques  such  as  those  given  in  Reference  14.  However,  F (<j>) 
is  not  known  initially,  but  it  is  assumed  that  we  have  a  reasonable  estimate  so  that  the  equations 
(5)  can  be  integrated  as  above.  Using  the  estimated  F (4>)  then,  after  integration,  a  “residual  normal 
velocity  given  by  the  value  of  the  left-hand  side  of  equation  (6)  at  £  =  0  will  be  obtained.  The 
problem  then  is  to  minimise  the  residual  at  £  =  0,  for  all  4>,  by  changing  the  function  F  (4>).  This 
minimisation  is  earned  out  by  the  iteration  procedure  described  in  Appendix  A. 

In  the  method  described  in  Appendix  A  it  is  necessary  for  the  conical  flow  solutions  required 
here  to  have  a  good  initial  estimate  of  F(<f>),  so  that  the  integration  will  not  diverge.  To  obtain  a 
good  estimate  of  the  shock  shape  F(4>)  for  conical  bodies,  the  following  procedure  is  adopted.  The 
flow  field  solution  is  first  found  for  a  circular  cone  at  zero  incidence  by  means  of  the  iteration  pro¬ 
cedure  described  in  Appendix  D.  A  very  small  perturbation  is  then  made  in  either  body  shape  or 
incidence  (for  example  an  elliptic  cone  of  small  eccentricity  or  a  circular  cone  at  very  small  incidence 
is  considered)  and  for  this  situation  a  solution  is  sought,  using  for  the  initial  estimate  of  F (<£)  the 
value  obtained  for  the  circular  cone  at  zero  incidence.  Having  obtained  the  solution  for  the  small 
perturbation,  a  much  bigger  perturbation  of  the  same  type  is  made  and  the  function  F(<)>)  extra¬ 
polated  so  that  a  good  initial  estimate  is  still  available.  For  example,  in  the  case  of  a  circular  cone 
at  incidence,  a  solution  is  first  found  for  an  incidence  a  where  a/0c  =  0.01,  then  F {<f>)  is  extrapolated 
linearly  to  a/dc  —  5a  where  5a  —  0.1  to  give  an  initial  estimate,  and  the  final  solution  at  that  in¬ 
cidence  is  found  by  iteration  as  mentioned  above.  After  this  computation  a  quadratic  extrapolation 
of  F (4>)  is  made  to  give  an  initial  estimate  at  a/6c  =  2  5a  and  so  on  to  higher  incidences. 

In  the  case  of  the  circular  cone  at  incidence,  the  function  F (4>)  is  represented  adequately 
by  a  F  ourier  series 

F  (<f>)  =  J2  Fi  cos  i  <f>  (8) 


since  F  (<£)  must  be  symmetrical  about  4>  =  0  and  <t>  = ir  (<j>  =  0  is  the  windward  plane  of  symmetry 
and  4>  =  ir  is  the  leeward  plane  of  symmetry).  This  form  has  the  advantage  that  m  can  be  kept 
small  ( =  1  or  2)  for  small  incidences,  thus  giving  fewer  “unknowns”  to  the  problem  and  economising 
on  the  iteration  process  to  determine  F  (4>).  The  value  of  m,  however,  must  be  increased  as  the 
incidence  is  increased.  Now  as  m  increases,  the  representation  (8)  of  the  function  F  (</>)  may  still 


be  good,  but  the  derivative  of  the  F ourier  series  may  not  be  a  good  representation  of 
j)F 

of  obtaining  by  differentiation  of  (8)  the  difference  formula  (A-l)  is  used. 

0(f> 


OF 

d(f> 


,  so  instead 
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A  difficulty  that  must  be  overcome  in  integrating  the  ordinary  differential  equations  that 
represent  the  set  (5)  .is  that  some  of  the  derivatives  are  not  defined  at  £.  =  0.  This  is  observed  by 
examining  the  equations  at  £  =  0  where  we  have 


„  ,  GG  .  W  A 

G  d  Uc  — —  P;  +  ^  us  =  0 

£r  d  vf  +  -  p5  +  q-  Vtf  =  q- 

.  I  G  Go  .  w  ,  po 

£rd  w«  ~—pGvt  +  Gw*  +  PG 


wv 

G 


.  ,  k  ..  ,  £o  pc2  ,  ,  pc-  .  w 

G  pc-  u5  +  G  Pc2  v5  H - w5  +  d'  P5  +  ^  w0  +  —  p0 

G  p  u5  +  G  p  vs  +  w5  +  d'  p5  +  -  wfl  +  —  p,i  =  —  - 


pc2  v 
r 


(a) 

(b) 

(c)  (9) 

(d) 

(e) 


W  5G"  £ 

where  d  =  v  —  G  u  —  — ,  d'  =  G  u  +  G  v  +  —  w,  and  the  subscripts  z,  r,  6,  £  and  $  refer 

•  (jt  d(p  r 

to  partial  differentiation  with  respect  to  these  variables. 

It  is  convenient  to  rewrite  (9d)  and  (9e)  in  terms  of 

a 


^G  U  +  £rV+y  wj 


(10) 


where  G>  G>  and  G  are  given  by  equation  (4)  with  z  =  1.  Equations  (9d)  and  (9e)  were  written  above 
without  the  limit  of  £  =  0,  thus  considering  the  derivative  (10)  in  connection  with  (9d)  and  (9e)  is 
permissible. 


Now  (10)  can  be  written 


w  /  F(5  —  G(1\ _ |o 

r  \  ¥  -  G  )  r2 


w  (F  —  G) 


G  +  £rvf  +  ~  w{  —  u 
and  (9d)  and  (9e)  as 

3°'  u  +  S* v  +_7  w)  +  u  +  G  +  G2  w  (F  “  G)] 


pc- 


k  J  ,  pc2  .  W  pc2  V 

+  £r  Cl  p£  +  -jrT  ~  Q 


(9d') 


[|  (*‘  u  +  G  +  7  w )  +  u  +  5  (f’  -  G ')  +  G  w  (F  G)] 


,  t  ,1  ,  PW0  W  £0  _  pv 

+  Gdp5  +.-Q-+  -G-  “  G 


(9e') 


on  applying  the  limit  £  =  0. 
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The  boundary  condition  at  the  surface  (eq.  (6)  )  implies  that  d  =  0  at  £  =  0  and  hence  it 
is  not  possible  to  determine  m,  v=,  w,,  or  p=  at  £  =  0  from  the  surface  equations  (9).  However,  if  we 
assume  that  du?,  dv=,  dw;,  dp,,  and  dp--  all  tend  to  zero  as  £  -»■  0,  which  appears  reasonable,  and.  also 
assume  that  derivatives,  of  quantities  with  respect  to  <f>  are  finite  on  the  surface  (which  is  reasonable 
since  the  surface  is  a  streamline),  then  we  can  at  least  say  that  p5  is  finite  from  (9a),  (9b),  or  (9c), 

and  that  ^  £,  u  +  £r  v  +  ~  w 

The  previous  analysis  shows  that  integration  of  the  differential  equations  right  up  to  £  =  0 
cannot  be  made  by  the  usual  techniques,  since  these  would  require  calculation  of  all  derivatives 
with  respect  to  £  at  £  =0.  However,  integration  can  be  made  to  a  value  <5£  close  to  the  surface  and 

t 

the  variable  £„  u  +  £r  v  +  ~  w  can  be  extrapolated  from  its  value  and  derivative  at  points  near 

to  the  surface  (at'  5£,  2  6£,  3  6£,  and  4  5£,'  for  example).  Extrapolation  of  this  quantity  is 
permitted,'  since  we  have  shown  its  derivative  is  at  least  finite  at  £  —  0.  Now  the  quantity 

£z  u  +  £r  v  +  -w  is  the  normal  velocity  V„  at  the  surface,  which  we  are  trying  to  make  zero 
1  -1  5  =  o  ° 

by  the  iteration  procedure.  Thus  the  problem  is  solvable  by  the  technique  mentioned. 

Having  completed  the  iteration  to  make  normal  velocity  Y  sufficiently  small,  a  final 
integration  of  the  set  (5)  is  made  to  a  value  of  £  fairly  close  to  the  surface  (approximately  0,003 ) 
and  then  pressure  is  obtained  at  the  surface  £  =  0  by  inward  extrapolation  from  0.003.  Again  this 
extrapolation  is  permitted  since  the  derivative  of  p  is  finite  at  the  surface.  The  values  of  Vu.  and  p 
thus  obtained  at  the  surface  allow  a  complete  calculation  of  the  other  flow  variables  at  the  surface,,  as 
follows. 


J  from  (9d')  or  (9e')  is  also  finite  at  £  =  0. 


Since  the  body  is  a  streamline  we  have 

^  )  =  constant  (11) 

P  /  c 

and  the  constant  is  equal  to  the  value  at  a  saddle  point  or  points  of  attachment,  on  the  surface 
which  for  the  circular  cone  at  incidence  is  at  the  windward  axis  <p  =  0  since  this  axis  also  forms  the 
same  streamline,  as  pointed  out  in  Reference  7. 

Also  in  the  flow  field,  and  hence  on  the  surface  for  isoenergetie  flow 

I  (u2  +  v2  +  w2)  =  constant  (12) 

7-1  p 

where  the  constant  is  determined  from  free  steam  conditions. 

A  further  condition  follows  by  eliminating  p5  from  (9a)  and  (9b)  (with  d  =  0) 

w  u„  +  w  G  v„  =  w2  G  (13) 

The  equations  (11),  (12),  (13),  together  with  the  boundary  condition  (6)  and  the- knowledge 
of  the  pressure,  give  five  equations  to  be  solved  for  u,  v,  w,  p,  and  p  at  the  surface.  A  method  of 
solving  these  equations  for  the  circular  cone  at  incidence  is  given  in  Appendix  B.  The  same  method 
can  be  used  for  other  conical  bodies,  provided  the  positions  of  the  saddle  points  of  attachment  on  the 
surface  are  first  determined. 
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It  should  be  noted  that  a  complete  surface  solution  may  not  always  be  necessary  since 
the  pressure  is  found  correctly  by  extrapolation,  and  it  is  probably  better  to  use  extrapolated  or 
near  the  surface  values  of  other  quantities  for,  say,  boundary-layer  calculations.  This  follows  from 
the  assumption  that  the  vortical  layer  (the  layer  near  the  surface  where  large  gradients  are  present) 
is  much  thinner  than  the  boundary  layer. 

4.0  A  SHORT  DISCUSSION  OF  THE  VORTICAL  LAYER 

It  has  been  pointed  out  (Ref.  7,  8,  and  9)  that  in  non-axisymmetric  conical  supersonic  flow, 
large  gradients  in  some  quantities  are  present  very  near  to  the  surface.  It  was  shown  in  Section  3.0 
by  examining  the  equations  of  motion  on  the  surface,  that  ,  the  derivatives  of  pressure  p  and  normal 

velocity  £z  u  T  £r  v  +  %  w  were  finite  at  the  surface,  while  derivatives  of  other  quantities  were 
indeterminate.  Since  only  the  normal  velocity  at  the  surface  was  required  for  iteration  and  the  pres¬ 
sure  p  was  required  to  solve  the  surface  equations  completely,  the  present  method  overcame  any 
difficulties  due  to  large  gradients  within  the  vortical  layer.  For  completeness,  however,  an  inte¬ 
gration  of  the  set  (5)  was  performed,  taking  very  small  increments  in  £  as  the  surface  was  approached.. 
Figure  2  shows  plots  of  the  velocity  components  u  and  w  and  the  density  profile  near  to  the  surface 
for  the  particular  case  of  a  circular  cone  with  Moo  =  7,  A  =  25,  a  =  10  at  <j>  =  90°. 

Also,  for  completeness,  Table  1  shows  a  .  comparison  between  surface  values  of  pressure 
obtained  by  integration  near  to  the  surface  and  those  obtained  in  the  usual  method  of  this  Report 
(i.e.  by  inward  extrapolation  from  £  «  0.003).  The  tabulated  results  are  for  the  case  Mco  =  1.797, 
A  =  12.5°,  a  =  7.5°;  this  example  is  chosen  since  there  is  a  fairly  noticeable  change  in  pressure 
near  to  the  surface  (for  the  above  case  of  Mco  =  7,  A  =  25,  «  =  10°  there  is  very  little  change  in 
pressure).  The  adequacy  of  the  usual  method  of  extrapolation  is  well  illustrated. 

5.0  COMPARISON  OF  THE  PRESENT  RESULTS  WITH  THOSE  OF  OTHER  THEORETICAL  METHODS  FOR 
THE  CIRCULAR  CONE  AT  INCIDENCE 

Figures  3,  4,  and  5  show  comparisons  of  the  present  theory  with  the  results  of  first  order 
theory  given  in  Reference  4.  Figure  4  also  shows  a  comparison  with  the  results  of  Moretti  (Ref.  13), 
while  Tables  2  to  9  give  comparisons  of  the  present  theory  with  that  due  to  Babenko  et  al.  (Ref.  11). 

It  can  be  seen  from  these  Tables  that  almost  exact  agreement  is  obtained  between  the 
present  results  and  those  of  Reference  11,  while  the  results  of  Reference  13  differ  by  a  larger  amount, 
probably  because  of  the  differences  in  step  size  used. 

The  step  sizes  used  by  Babenko  were  5£  =  0.05  and  S<f>  =  11.25°,  by  Moretti  S£  =  0.167 
and  8<j>  =  18°,  and  with  the  present  method  5£  =  0,1  and  o<p  —  22^°  for  most  of  the  cases  in 
which  comparisons  were  made  in  this  Report.  Good  accuracy  is  obtained  with  the  present  method 
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even  though  larger  step  sizes  are  used.  This  is  to  be  expected  since  the  errors  are  0  (A/i )  and  0  (5 £  ) 
in  the  present  scheme. 

The  first  order  results  agree,  as  expected,  for  small  relative  incidences  only,  and  the  inade¬ 
quacy  of  the  first  order  theory  is  noted  even  for  relative  incidences  of  about  0.5.  Results  from 
References  11  and  13  are  available  cnly  for  relative  incidences  a  /A  up  to  about  0.8.  At  present, 
accurate  tabulated  results  for  higher  relative  incidences  are  not  available,  so  that  comparisons  cannot 
be  made  with  other  theoretical  methods  for  relative  incidences  higher  than  0.8.  However,  a  compar¬ 
ison  with  experimental  results  is  made  in  the  next  Section,  which  extends  to  higher  relative  incidences. 


6.0  COMPARISON  OF  PRESENT  RESULTS  WITH  EXPERIMENT  FOR  THE  CIRCULAR  CONE 

Experiments  to  measure  surface  pressures  on  12.5°  and  5°  circular  cones  have  been  made 
at  NAE  under  high  Reynolds’  number  conditions  (Ref.  15  and  16).  The  experiments  were  con¬ 
ducted  over  a  wide  range  of  Mach  numbers  from  1.8  to  4.25  and  with  relative  incidences  from  zero 
to  about  2.5.  Figures  6  —  9  show  a  comparison  of  surface  pressure  ratio  p/pco  as  calculated  by  the 
present  method  with  the  experimental  data. 

In  most  cases  the  greatest  difference  between  the  theory  and  experiment  occurs  in  the 
region  <j>  =  0°  to  30°.  It  is  known  that  the  high  surface  shear  stress  in  this  region  leads  to  appreciable 
positive  hole  errors  in  the  measured  static  pressures  (Ref.  17),  but  the  exact  magnitude  is  difficult 
to  estimate  in  each  case. 

It  can  be  seen  from  the  Figures  that  excellent  agreement  is  obtained  in  the  region  from 
<j>  =  30°  to  180°,  even  at  incidences  where  the  singularity  has  left  the  surface,  e.g.  =  1.8, 
0O  =  12.5°,  a  —  17.5°,  and  also  where  locally  conically  supersonic  conditions  are  just  becoming  present 
near  the  surface,  e.g.  =  4.25,  9C  =  12.5°,  a  =  12.5°.  The  present  method  is  limited  to  inciden¬ 
ces  below  some  critical  value_at  which  a  region  of  the  ffow  becomes  locally  conically  supersonic, 
However/  it  can  be  slen  that  surface  pressures  extrapolated  from  incidences  below  this  critical 
incidence  (by  quadratic  extrapolation)  still  give  good  agreement  with  the  experiment,  at  least  in 
the  region  of  the  flow  where  separation  does  not  have  large  effect  (see,  for  example,  Fig.  10). 

7.0  THE  ELLIPTIC  CONE  IN  SUPERSONIC  FLOW 

Most  of  the  discussions  in  this  Report  have  been  concerned  with  circular  cones  at  incidence. 
However,  the  method  of  solution  is  applicable,  at  least  in  principle,  to  any  conical  body  with  coni- 
cally_subsQliic_flow.  As  a  further  example  of  the  method  some  calculations  were  made  for  elliptic 
cones  at  zero  incidence,  and  at  incidence  but  without  yaw,  again  starting  from  the  a  =  0  circular 
cone  solution.  Some  of  these  results  (the  surface  pressures)  are  compared  with  another  theoretical 
method  (Ref.  18)  and  with  experiments  (Ref.  19)  in  Figures  11  and  12.  It  can  be  seen  from  these 
Figures  that  the  present  method' gives  much  better  agreement  with  the  experiments,  particularly 
for  the  cases  at  incidence,  than  the  linearised  characteristics  solutions  of  Martellucci  (Ref.  18). 

The  method  has  been  found  to  be  less  efficient  when  the  ratio  of  major  to  minor  axes, 
a/b,  becomes  so  high  that  large  gradients  in  quantities  occur  in  the  flow  field  near  the  “leading 
edge”.  Very  small  increments  5<j>  are  needed  for  the  difference  scheme'  (A-l).  A  reasonable  limit  for 
the  ratio  a/b  seems  to  be  within  the  range  2  to  3,  (see  computer  times,  App.  E.4.2). 

8.0  A  FURTHER  EXAMPLE  OF  A  CONICAL  BODY  IN  SUPERSONIC  FLOW 

As  a  further  example  to  illustrate  the  use  of  the  present  method,  computations  were  made 
for  the  conical  body  given  by 

,  r  =  G0  +  G2  cos  2<j>  +  G4  cos  4 <f> 

at  a  unit  distance  z  =  1  from  the  body  tip,  where  G0,  G2,  and  G4  are  constants  (0.2679,  —0.01,  and 
0.02  respectively,  in  the  example  for  which  results  are  given).  The  computation  was  started  from  a 
circular  cone  at  zero  incidence  (G2  =  G4  =0)  and  G2  and  G4  were  given  small  values  (e2,  e4,  say) 
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and  a  solution  obtained,  G2  and  G.i  were  then  increased  proportionally  to  e2  and  e4  and  an  initial  es¬ 
timate  of  the  shock  shape  was  made  by  extrapolation  from  the  solution  with  G2  =  e2  and  G4  =  a, 
and  a  solution  was  again  obtained  by  iteration.  G2  and  G4  were  increased  again  proportionally  to 
62  and  64  and  the  process  continued. 


Figure  13  shows  the  shock  shape  and  the  surface  pressure  distribution  corresponding  to 
Go  =  0.2679,  G2  =  —0.01,  G4  =  0.02  at  zero  incidence.  Figure  14  shows  pressure  distributions  and 
also  indicates  the  circumferential  angles  at  which  nodal  and  saddle  type  singular  points  on  the 
surface  are  present  for  the  same  body  at  several  incidences. 


This  example  illustrates  well  the  singularity  behaviour  to  be  expected  for  this  type  of 
indented  body.  It  can  be  seen  that  at  zero  incidence  there  are  two  nodal  and  three  saddle  type 
singular  points  on  the  surface  in  the  range  0  <  <$>  <  180.  As  the  incidence  is  increased  the  nodal  and 
one  of  the  saddle  type  singularities  on  the  windward  side  combine  to  “cancel”  each  other  so  that  the 
surface  streamlines  now  converge  only  near  to  the  leeward  generator.  As  the  incidence  is  further 
increased  the  nodal  type  singularity  near  to  the  leeward  generator  moves  around  to  eventually  lie 
on  the  leeward  generator  where  previously  there  was  a  saddle  type  singular  point,.  The  positions 
of  the  singularities  are  easily  determined  from  the  direction  of  the  streamlines  as  integration  is  made 
into  the  surface,  and  it  is  also  known  that  they  can  be  located  at  points  on  the  surface  only  where  the 
circumferential  gradient  of  pressure  is  zero. 


9.0  CONCLUSIONS 

A  method,  has  been  presented  that  shows  that  the  elliptic  partial  differential  equations 
defining  the  conical  inviscid  flow  between  a  conical  body  and  its  conical  shock  wave  are  solved  very 
efficiently  by  a  numerical  approach. 

A  transformation  is  used  that  has  the  effect  of  including  the  body  shape  in  the  coefficients 
of  the  differential  equations  and  of  the  boundary  conditions,  thus  making  the  computer  program 
suitable,  in  principle,  for  any  conical  situation  simply  by  writing  a  few  program  statements  to  define 
the  body  shape. 

Solutions  have  been  obtained  by  making  successive  perturbations  to  a  circular  cone  at 
zero  incidence  for  which  the  flow  field  solution  is  readily  available.  Perturbations  of  incidence  were 
made  to  obtain  results  for  circular  cones  at  incidence  and  have  given  good  results  even  at  relative 
incidences  higher  than  any  other  known  fully  numerical  method.  To  generate  solutions  for  an 
elliptic  cone,  perturbations  were  made  first  in  body  shape  and  then  in  incidence,  if  solutions  at 
incidence  were  required.  Similarly,  solutions  can  be  obtained  for  other  conical  body  shapes,  such  as 
those  obtained  for  the  body  given  by  a  fourth  order  even  cosine  Fourier  series. 


It  was  shown  that  excellent  agreement  with  other  theoretical  methods '(in  particular  that 
due  to  Babenko  et  al.  (Ref.  11)  for  the  circular  cone)  and  also  with  experiment,  was  obtained.  At  the 
same  time,  solutions  were  generated  at  about  30  to  50  times  faster  with  the  present  method  (typical 
computer  times  are  quoted  in  App.  E.4.0)  than  with  other  fully  numerical,  accurate  techniques. 
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At  present  the  method  appears  to  be  limited  to  cases  that  do  not  have  regions  of  conically 
supersonic  flow,  since,  if  such  regions  do  exist,  the  defining  equations  in  these  regions  become  hyper¬ 
bolic  while  they  remain  elliptic  in  the  other  regions.  Two  other  conditions  also  present  a  limit  to 
the  applicability  of  the  method.  These  occur  when  an  entropy  singularity  moves  sufficiently  far 
away  from  the  surface  that  it  lies  very  near  to  the  first  exterior  set  of  mesh  points,  and  also  when  the 
shock  wave  approaches  very  close  to  the  Mach  cone  from  the  apex,  which  occurs  for  very  slender 
conical  flows  at  low  supersonic  Mach  numbers. 
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TABLE  1 

p  NEAR  TO  THE  SURFACE  FOR  THE  CIRCULAR  CONE 

Mro  =  1.797,  dc  =  12.5°,  «  =  7.5° 


£ 

0 

22.5 

45 

67.5 

90 

112.5 

135 

157.5 

180 

0.1 

0.8473 

0.8322 

0.7919 

0.7382 

0.6857 

0.6463 

0.6230 

0.6124 

0.6096 

0.05 

0.8522 

0.8361 

0.7930 

0.7359 

- 

0.6809 

0.6412 

0.6197 

0.6116 

0.6099 

0.025 

0.8537 

0.8370 

0.7924 

0.7333 

0.6769 

0.6371 

0.6169 

0.6107 

0.6099 

0.0125 

0.8541 

0.8371 

0.7916 

0.7315 

0.6743 

0.6346 

0.6152 

0.6101 

0.6099 

0.00625 

0.8542 

0.8370 

0.7911 

0.7304 

0.6729 

0.6331 

0.6142 

0.6098 

0.6099 

0.003125 

0.8543 

0.8369 

0.7908 

0.7299 

0.6721 

0.6324 

0.6137 

0.6096 

0.6099 

1.53xl0~6 

0.8543 

0.8369 

0.7905 

0.7293 

0.6713 

0.6316 

0.6132 

0.6095 

0.6099 

Value 

extrapolated 

from 

$  =  0.003125 

0.8543 

0.8369 

0.7905 

0.7293 

0.6713 

0.6316 

0.6132 

0.6095 

0.6099 

f 


TABLE  2 


COMPARISONS  OF  SURFACE  PRESSURE  AND  SHOCK  SHAPE  BETWEEN  PRESENT 
THEORY  AND  THE  THEORY  OF  BABENKO  ET  AL.  (REF.  11) 

Mm  =  5 ,  0  —  10,  a  =  7.5 


0 

22.5 

45 

67.5 

.6170 

0.5872 

0.5081 

0.4046 

li  •' 

.6172 

0.5874 

0.5082 

0.4047 

.2488 

0.2507 

0.2567 

0.2674 

.2489 

0.2508 

0.2568 

' 

0.2674 

90 

112.5 

135 

157.5 

0.3044 

0.2284 

0.1894 

0.1822 

0.3045 

0.2284 

0.1895 

0.1821 

0.2832 

0.3041 

0.3281 

0.3499 

0.2831 

0.3040 

0.3281 

0.3495 

TABLE  3 


Subscripts:  J  Values  obtained  by  present  method 


B  Values  obtained  by  Babenko  et  al.  (Ref.  11) 


TABLE  4 


COMPARISONS  OF  SURFACE  PRESSURE  AND  SHOCK  SHAPE  BETWEEN  PRESENT 
THEORY  AND  THE  THEORY  OF  BABENKO  ET  AL.  (REF.  11) 

=  2,  ec  =  15,  «  =  10 


4> 

22.5 

"  1 

45 

67.5 

90 

112.5 

135 

157.5  j 

■ 

Pj 

1.0143 

0.9839 

■ 

0.9018 

0.7931 

0.6884 

0.6153 

0.5833 

0.5783 

. 

Pb 

1.0153 

0.9852 

0.9039 

0,7957 

0.6914 

0.6178 

0.5852 

0.5805 

Fj 

0.5624 

0.5691 

0.5896 

0.6244 

0.6737 

0.7351 

0.8010 

0.8550 

Fb 

0.5636 

0.5703 

0.5909 

0.6259 

0.6755 

0.7369 

0.8037 

0.8574 

_ i 

TABLE  5 

Mm  =  5,  0O  =  15,  a  =  10 


0 

22.5 

45 

67.5 

90 

112.5 

135 

157.5 

pj 

1.0141 

0.8667 

0.6750 

] 

0.4897 

Pb 

1.0144 

0.8668 

0.6748 

0.3467 

Fj 

0.3459 

0.3476 

0.3530 

0.3623 

0.3942 

Fb 

0.3460 

0.3477 

0.3530 

0.3624 

0.3763 

1 

0.3943 

_ 1 

Subscripts;  J  Values  obtained  by  present  method 


B  Values  obtained  by  Babenko  et  al.  (Ref.  11) 


TABLE  6 


COMPARISONS  OF  SURFACE  PRESSURE  AND  SHOCK  SHAPE  BETWEEN  PRESENT 
THEORY  AND  THE  THEORY  OF  BABENKO  ET  AL.  (REF.  11) 

M#  =  7,  0e  =  15,  a  =  10 


0 

22.5 

45 

67.5 

90 

112.5 

135 

157.5 

180 

1.0795 

1.0178 

0.8542 

0.6435 

0.4426 

| 

0.2901 

0.1972 

0.1615 

0.1562 

.0798 

1.0179  . 

0.8544 

0.6433 

0.4426 

0.2899 

0.1974 

0.1615 

0.1560 

).3260 

0.3271 

0.3363 

0.3561 

0.3681 

0.3755 

0.3772 

J.3261 

0.3272 

0.3363 

0.3449 

0.3561 

0.3679 

0.3753 

0,3776 

Subscripts:  J  Values  obtained  by  present  method 

1) 


TABLE  7 


COMPARISONS  OF  SURFACE  VALUES  AND  SHOCK  SHAPE  BETWEEN  PRESENT 
THEORY  AND  THE  THEORY  OF  BABENKO  ET  AL.  (REF.  11) 

Mo,  =  5,  dc  =  25,  a  =  20 


4> 

22.5 

45 

67.5 

90 

112.5 

135 

157.5 

180 

Uj 

1.3167 

1.3571 

1.4205 

1.5026 

1.5956 

1.6906 

1.7741 

1.8123 

uB 

1.3026 

1.3165 

1.3572 

1.4217 

1.5048 

1.5989 

1.6936 

1.7711 

1.8129 

Vj 

0.6328 

0.6624 

0.8273 

0.8451 

VB 

0.6139 

0.7456 

0.8259 

0.8454 

Wj 

0.1785 

0.4862 

0.6358 

0.4319 

0 

wB 

0 

0.1792 

0.3446 

0.4832 

0.5818 

0.6285 

0.4463 

0 

Pj 

2.6838 

2.5058 

1.4619 

0.5447 

0.2522 

PB 

2.6842 

2.5062 

2.0423 

1.4596 

0.9282 

0.5434 

0.2508 

Pj 

4.7758 

4.5474 

3.9284 

2.2362 

1.5288 

0.8580 

0.8819 

Pb 

4.7759 

4.5475 

3.9289 

2.2368 

0.8603 

0.8785 

Fj 

0.5943 

0.6168 

0.6657 

0.6920 

fb 

0.5947 

0.6173 

0.6388 

0.6665 

0.6949 

0.6917 

Subscripts:  J  Values  obtained  by  present  method 


B  Values  obtained  by  Babenko  et'al.  (Ref.  11) 


SHOCK  WAVE 


CO-ORDINATE  SYSn 


FIG.  3  ••  CIRCULAR  CONE 

PLOTS  OF  FIRST  ORDER  AND  PRESENT  SURFACE  PRESSURE  RESULTS 

7  ec=25° 


FIG.  4  :  CIRCULAR  CONE 

PLOTS  OF  FIRST  ORDER,  MORETTl'S,  AND  PRESENT  SURFACE 

PRESSURE  RESULT S 


FIG.  6 '■  CIRCULAR  CONE 

COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  PRESENT  THEORY  AND  EXPERIMENT 

Moo=  4-25,  6c=5°,  a  =5° 


FIG.  7  :  CIRCULAR  CONE 

COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  PRESENT  THEORY 

AND  EXPERIMENT 


FIG.  9  :  CIRCULAR  CONE 

COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  PRESENT  THEORY 

AND  EXPERIMENT 

M  =4  25,  ec  =  l2  5°.  a  =  8-24°  AND  12-5° 


A  EXPERIMENT— ZAKK AY  AND  VISICH  (REF  19) 

-  THEORY— .MARTELLUCCI  (REF  18) 

X  PRESENT  THEORY 


FIG.  II'.  ELLIPTIC  CONE 

COMPARISON  OF  SURFACE  PRESSURE  BETWEEN  PRESENT  THEORY, 
THEORY  OF  MARTELLUCCI,  AND  EXPERIMENT 

fa b  =  0-3017  ,  a  =  0° 


FIG.  13:  BODY  GIVEN  BY  r  =  G0+  G2  COS  2  4>  +  G4C0S  4  <j>  (G0  =  0-2679,  G2=  -O  OI,  G4=0  02) 


SHOCK  SHAPE  AND  PRESSURE  DISTRIBUTION  FOR  M=  5  ,  a  =  0' 


SOLUTION  STARTED  FROM  CIRCULAR  CONE  r  =  Gc 
ALTERED  BY  CHANGING  G„  ,  G„ 


AND  BODY 


WINDWARD 


LEEWARD 


FIG.  1 4 '•  BODY  GIVEN  BY  r  =  G0+G2COS  2<£  +  G4COS4<£  (G0=  0-2679,  G2= -001,  G4=  0  02  AT 

UNIT  DISTANCE  FROM  APEX) 

M  =5,  SURFACE  PRESSURE  DISTRIBUTIONS  FOR  VARIOUS  INCIDENCES 


FIG.  15-  CIRCULAR  CONE 

PLOTS  OF.D  =  .j-  (£u+£v)  AND  ALONG  LEEWARD  PLANE  OF  SYMMETRY  (<£  =  7 r) 

ENTROPY  SINGULARITY  IS  SITUATED  WHERE  ^(£zu+£rv)  =  0 

M„=  I-797-  6  =12-5° 


FIG.  16:  CIRCULAR  CONE 

PLOT  OFp^+2Jov2AT  THE  LEEWARD  GENERATOR  </>■=  it  (USING  ISENTROPIC 
SURFACE  VALUES)  AGAINST  RELATIVE  INCIDENCE  -f  FOR  THE  CASE 

.  .  M  =  1-797  8  =12  5° 
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FIG.  17  :  ELLIPTIC  CONE 

PRESSURE  DISTRIBUTION  AROUND  ELLIPTIC  CONES  AT  ZERO 

INCIDENCE 

Meo=  10  a  =  0-7  ' 
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APPENDIX  A 

NUMERICAL  PROCEDURE  FOR  SOLVING  ELLIPTIC  PARTIAL  DIFFERENTIAL  EQUATIONS 


A.  1.0  NOTATION 

(x,  y)  co-ordinate  system 

U  a  function  of  (x,  y) 

x0,  xx  lower  and  upper  limits  on  x 

y<b  yi  lower  and  upper  limits  on  y 

F  (y)  a  function  of  y  at  x  =  x0 

e  (y)  a  function  of  y  at  x  =  Xj 

Xi 

x0 

Pi  a  value  of  p  on  the  ith  line 


ek  value  of  e  (y)  at  the  kth  line  at  x  = 

F  i  value  of  F  (y)  at  the  ith  line  at  x  = 

U  i  a  value  of  U  on  the  ith  line 


boundary  condition  at  x  =  Xi  is  that  f 


A. 2.0  METHOD 


Suppose  '(x,  y)  are  independent  variables  and  U  is  any  function  of  (x,  y)  defined  by  partial 
differential  equations  within  a  region  x0  <  x  <  Xi  and  y„  <  y  <  yi.  Boundary  conditions  for  U  or 
its  normal  derivatives  are  given  on  the  bounds  of  the  above  region. 

By  estimating  some  unknown  function  (or  functions)  at  one  of  the  boundaries  (say,  for 
dU 

example,  —  =  F  (y)  at  x  =  x0,  y0  <  y  <  yi),  and  by  replacing  the  derivatives  in  the  y  direction  by 
differences,  thus  making  the  partial  differential  equations  into  ordinary  differential  equations,  the 
equations  can  be  integrated  from  x0  to  Xi  in  a  way  similar  to  that  for  integrating  parabolic  partial 

differential  equations.  At  x  =  x,  there  are  given  boundary  conditions,  of  the  form  f  ^U,  ^  j  =  0, 
to  be  satisfied,  but  the  preceding  integration,  assuming  the  estimate  at  x  =  x0  is  not  correct,  would 

=  e  (y),  say.  To  solve  the  elliptic  problem  completely  the  estimated 
function  at  x  =  xp  must  be  improved  until  e  (y)  is  sufficiently  small. 


give  a  residual  of  f  ^U, 


dU 

dX 
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dll 


x0)  .the  region  is  first 
divided  into  strips  of  width  8y  and  first  or  second  derivatives  replaced  by  differences 


To  carry  out  the  above  integration  (using  an  estimate  of  at  x 


=  U  (x,-  y  +  Sy)  -  U  (x,  y  -  8y)  +  Q  (gy2) 


u  (x,  y  +  8y)  -  2U  (x,  y)  +  U  (x,  y-Sy) 

sy2 .... 


+  o  (ay*) 


or  alternatively  by  more  accurate  formulae  such  as 


\  4  riJ  (x,  y  +  Sy)  -  U  (x,  y  -  6y)' 

1  fU  (x,  y  +  28y)  -  .U  (X,  y  -  25y)l 

/  x,y  3  L  2  Sy 

j 

1  >1 
1*0 

|co 

i 

T  0  (ay4) 


(A-l) 


for  a  first  derivative.  The  partial  differential  equations  thus  become  a  set  of  coupled  ordinary 
differential  equations  in  one  independent  variable  x  with  differential  equations  for  eaeh  of  the 
dividing  lines  in  the  region  y0  :<  y  ^  y,  and  with  the  differential  coefficient  at  any  line  depending 
on  variables  to  both  sides  of  that  line.  These  resulting  ordinary  differential  equations  can  then  be 
integrated  by  standard  techniques,  one  of  the  most  efficient  being  the  Hamming  predictor  modifier 
corrector  method  (Ref.  14).  This  method  was  used  in  calculations  presented  in  this  Report;  the 
starting  procedure  employed  was  the  Runge  Kutta  method  given,  for  example,  in  Reference  20. 

Once  an  integration  has  been  made  from  one  boundary  to  the  other  (x0  to  xx  in  the.  case 
above)  the  residual  function  e(y)  is  known  for  the  given  estimate  F(y).  To  improve  the  estimate 
F(y)  so  that  |e(y)|  is  made  smaller,  the  following  method  is  used. 


The  function  F(y)  can  be  defined  by  its  values  at  F(y0  +  j<5y)  (j  =  0,  1  .  .  .  n;  y0  +  n5y  = 
yx),  and  similarly  e(y)  is  represented  by  its  values  at  e(y0  +  joy).  Therefore  determining  a  procedure 

n 

to  minimise  ]e(y)|  is  equivalent  to  finding  a  procedure  to  minimise  e-  (y0  +  koy)  with  respect  to 

.'  _  k  =  °  ' 

F  (yo  +  j6y),  j  =  0,  1  ....  n.  Many  methods  exist  for  minimisation;  one  of  the  best  for  minimising 
a  sum  of  squares  is  presented  by  Powell  (Ref.  21).  This  method  is  similar  to  the  generalized  least 
squares  technique  given  by  the  iterative  process 

{k?o  ST,  Wi }  3Fi  =  ”  k?0  Ck  (i  =  0,  1  ...  m)  (A’2) 

where,  in  the  above  case,  the  number  of  unknowns  m  is  equal  to  n.  ek  =  e  (y0  +  k5y.)  and 


£ 
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Fj  ss  F  (y0  +  iSy)  and  5Fj  is  the  improvement  to  be  made  to  F,-  so  that  V  t  is  minimised. 

k  =  o  k 

In  the  generalised  least  squares  technique  “steps”  defined  by  (A-2)  are  made  until  either  e* 
or  SFj/Fj  is  sufficiently  small.  Each  step  requires  calculating  3ek/3Fj  by  differences 


d*k  ek  (Fp,  Fx  ,  ,  ,  ,  Fj  -f-  AFj,  .  .  .  Fm)  —  eic  (Fp,  Fi  .  .  .  Fj  .  .  ,  Fm)  (A-3) 

3Fi  “  AF, 

for  i  =  0,  1,  2  .  .  .  m, 


where  ek  is  considered  as  a  function  of  F0,  Fi.  .  .  Fm  since,  for  given  values  of  F0,  Fi  .  .  .  Fm,  the 
corresponding  values  of  ek  (k  =  0,  1  .  .  .  n)  can  be  found  by  integration  as  described  previously. 
AFi  is  a  small  increment  in  F„  say  10  6  Fi  if  F,  A  0.  Now  in  Powell’s  method  only  the  first  “step” 
requires  the  use  of  (A-3),  and  after  the  first  step  the  partial  derivatives  can  be  calculated  from  values 
of  ek  (k  =  0,  1  .  .  .  n)  already  obtained  on  the  previous  step.  Thus  Powell’s  method  is  more  efficient 
than  the  generalised  least  squares  technique  and  it  is  also  claimed  to  ensure  convergence. 

It  should  be  noted  that  it  is  not  necessary  to  define  the"  function  F(y)  by  its  values  at 
yo  +  joy  (j  =  0,  1  .  .  .  n).  In  fact  it  is  more  economical  to  define  F(y)  by  as  few  “unknowns”  as 
possible  in  order  to  reduce  the  computation  required  for  the  first  step  given  by  (A-3).  For  instance, 
the  first  few  terms  of  a  Fourier  series  expansion  could  be  used  if  it  is  known  that  F(y)  can  be  ade¬ 
quately  represented  in  this  way,  which  is  the  case  for  a  circular  cone  at  incidence  and  for  the  other 
conical  bodies  discussed  in  this  Report. 

A  simple  example  illustrates  the  above  method.  To  solve 


d2U  cPU  1 
dx2  +  dy2  +  1 


(A-4) 


given  U  =  0  on  x=  0,  1  and  on  y  =  0,1.  From  symmetry  it  is  sufficient  to  solve  the  problem  for 
the  region  0<x<J,  0<y<i  with  boundary  conditions 


U  =  0  on  x  =  0  and  on  y  =  0 


3U 

dX 


=  0  on  x  = 


|  and  ^  =  O'  on  y 
2  dy  y 


1 

2 
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Divide  the  region  0  <  x  <  2  ,  0  <  y  <  ^  into  strips  as  shown,  and  number  the  lines  forming 
the  strips  0,  1,  2,  ....  n 


c)^U 

Replace  —  2  at  line  i  by  the  difference  approximation 


Let  Pi 


-(£).• 


Ui+i  ~  2Uj  T  Ui-i 

ay2 

then  equation  (A-4)  can  be  written 


dpi  1  Ui+i  —  2U;  +  Ui_i  ,  1  A 
d*  +  - —  +1  =  0 


also 


dl JA 
dx 


Pi 


(A-5) 


Now  the  set  (A-5)  represents  a  set  of  ordinary  differential  equations.  The  boundary  condition 
dl] 

—  =  0  on  y  =  0.5  can  be  satisfied  by  the  difference  approximation  Un  + 1  =  Un  - 1.  The 

oy 

equations  (A-5)  are  subject  to  boundary  conditions  at  x  =  0 

Ui  =  0 


Pi  =  F(y,) 

On  estimating  F(y),  the  set  (A-5)  can  be  integrated  to  x  =  0.5  giving  residual  errors  e(ioy)  =  p, 
at  x  =  0.5  (we  require  p,  =  =  0  at  x  =  0.5).  One  can  then  make  small  displacements  to  F(y) 

as  in  (A-3)  and  use  (A-2)  to  improve  the  estimate  of  F(y). 

An  advantage  of  the  above  procedure  is  that  if  the  differential  equations  and  boundary 
conditions  are  linear  as  in  the  above  example,  then  the  correct  solution  will  be  obtained  after  only 
one  step  given  by  (A-2),  since  in  this  case  the  ek  are  linear  functions  of  F|.  A  disadvantage  of  such 
an  approach  is  that  reasonably  good  initial  estimates  of  F j  may  be  necessary  in  order  that  the  inte¬ 
gration  will  not  diverge. 
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APPENDIX  B 

A  METHOD  FOR  SOLVING  THE  SURFACE  EQUATIONS  FOR  THE  CIRCULAR 

CONE  AT  INCIDENCE 

The  surface  pressure  is  given  by  the  method  of  Section  3.0,  but  it  is  necessary  to  calculate 
u,  v,  w,  and  p  from  equations  (6),  (11),  (12),  and  (13).  Obviously  the  density  follows  directly  from 

equation  (11)  where  the  constant  in  that  equation  is  the  value  of  ^  at  <p  =  0. 

At  the  windward  and  leeward  axes  (<f>  =  0,  ir),  w  =  0  from  the  symmetry  of  the  problem, 
hence  at  these  two  points  u  and  v  can  be  calculated  from  (6)  and  (12) .  Equation  (13)  u(1  +  Gv(,  =  wG 
is  replaced  by  differences  as  in  equation  (A-l).  These  difference  equations  are  written  at  every 
point  where  p  is  given  (except  at  0,  v)  and,  together  with  (6)  and  (12)  written  at  every  such  point, 
give  a  set  of  non-linear  algebraic  equations.  These  equations  can  then  be  solved  by  usual  techniques. 
The  method  employed  here  was  that  of  Reference  21,  since  this  method  is  efficient  and  was  already 
being  used  to  solve  the  main  part  of  the  problem. 
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APPENDIX  C 

THE  POSITION  OF  NODAL  AND  SADDLE  TYPE  SINGULAR  POINTS 


It  has  been  noted  in  Reference  7  that  singular  points  will  always  exist  in  conical  flow 
without  axial  symmetry.  We  define  a  nodal  type  singular  point  as  a  point  where  surface  streamlines 
converge  and  a  saddle  type  singular  point  as  one  where  the  surface  streamlines  diverge.  At  a  nodal 
type  singular  point,  discontinuities  in  entropy  will  exist  giving  a  discontinuity  in  some  other 
quantities,  while  at  the  saddle  type  singular  point  no  discontinuities  will  exist.  In  this  Section  we 
will  consider  the  position  of  the  singular  points  for  a  circular  cone  at  incidence  and  also  for  an 
elliptic  cone  at  zero  incidence. 

The  position  of  a  singularity  is  important  since,  if  it  should  move  off  the  surface,  it  is  not 
possible  to  integrate  the  equations  (5)  past  this  singularity.  On  the  other  hand,  a  knowledge  of  their 
position  on  the  surface  is  also  important,  since  this  would  affect  the  solution  of  the  surface  equations 
(6),  (11),  (12),  (13).  .. 

We  know  that  a  singular  point  is  defined  by  the  two  conditions 

w  =  0  (C-la) 


and  the  main  diagonal  of  matrix  B  should  be  zero,  that  is 

£zu  +  £r  v  +  w  =  0 

since  equations  (9d)  and  (9e)  can  be  combined  to  give  (  £ru  +  £r  v  +  —  w)^  +  —  ^  = 

\  r  /of  v  otp 

S  =  log  -P  and,  hence,  when  conditions  (C-l)  are  satisfied  S  is  indeterminate. 


(C-lb) 
0  where 


Now  (C-lb)  is  satisfied  at  the  surface  as  a  boundary  condition,  and  because  of  the  symmetry 
about  </i  =  0  and  x  in  the  case  of  a  circular  or  elliptic  cone,  we  also  know  that  w  =  0  at  4> '  —  0  or  x 
(as  well  as  at  <f>  =  x/2  for  the  elliptic  cone  at  zero  incidence),  Thus  a  singular  point  is  always  present 
at  a  point  on  the  body  and  in  a  plane  of  symmetry  for  the  conditions  being  considered. 


C.1.0  THE  CIRCULAR  CONE  AT  INCIDENCE  ' 

It  is  shown  in  Reference  9  that  for  the  circular  cone  at  incidence  a  saddle  type  singular 
point  is  present  at  the  windward  generator  4>  =  0  and  that  a  nodal  type  singular  point  is  present 
at  the  leeward  generator  4>  —  x  until  the  surface  pressure  becomes  adverse.  It  is  then  shown  that 
this  nodal  singularity  may  move  away  from  the  leeward  generator,  to  follow  the  point  of  minimum 
circumferential  pressure,  or  that  it  may  stay  at  the  leeward  generator  until  a  certain  critical  incidence 
is  reached,  when  it  would  lift  off  the  surface.  The  following  analysis  shows  that  the  latter  situation  is, 
in  fact,  what  happens  to  this  nodal  singularity. 

To  determine  whether  the  nodal  singularity  can  move  along  the  leeward  plane  of  symmetry 
and  off  the  surface,  we  consider  the  quantity  [£.,.u  +  £r  v]^  _  1;  0  _  v  which,  as  the  normal  compo¬ 
nent  of  the  velocity  at  the  shock,  must  be  negative  in  the  sign  convention  being  used.  We  also  consider 


U(e-u  +  £'v)].,0 


0  “  7T 
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and  show  that  this  derivative  becomes  positive- under  the  conditions  given  in  Reference  9  for  the 
singularity  to  move  away  from  £  =  0,  4>  =  x.  Thus,  since  £2u  +  £rv  =  0  at  £  =  0  (the  normal 
velocity),  it  follows  that  £zu  +  £rv  will  pass  through  another  zero  along  the  leeward  plane  of  sym¬ 
metry  as  well  as  being  zero  on  the  surface. 

It  will  also  be  shown  that  the  nodal  singularity  will  not  move  along  the  surface  away  from 
the  leeward  generator.  Thus,  tire  nodal  singularity  will  move  into  the  flow  field  along  the  leeward 
plane  of  symmetry. 

C.1.1  To  Show  That  ~  {  £ZU  +  £rV  )  Becomes  Positive  At  A  Certain  Critical  Incidence 

\  /  -I  i  =  o,  V  =  X 

Consider  equation  (9c),  which  applies  at  the  surface  £  =  0 

£r  G0  w  p0  _  wv  ■ 

“^G~P5+  GwV+pG-  ~G_ 

since  d  =  0  is  a  surface  boundary  condition.  Now  differentiate  this  equation  with  respect  to  <f> 
and  consider  the  result  at  a  plane  of  symmetry  (G0  =  w  =  p0  =  0)  giving 

-  ;--G,<0  pE  +  w02  +  +  w0v  =  0 

p  p 


-  i  [  -  V  ±  J  *  -  4  (  f  -  T)  ] 


where  T  =  £r  G0(i  pE/p  and  must  be  zero,  since  from  (9b)  it  follows  that  if  w  =0  then  p5  =  0. 
For  the  circular  cone,  in  the  leeward  or  windward  plane  of  symmetry  the  positive  sign  obviously 
applies  at  incidence  a  =  0,  since  we  then  know  that  p!l(,  =  w(i  =  0  and  it  seems  from  the  results 
obtained,  using  the  method  of  this  Report  or  from  those  of  Reference  11,  that  the  -f-  sign  always 
applies  on  the  windward  side.  The  +  sign  applies  also  on  the  leeward  side  until  an  incidence  is 
reached  where  v2  =  4  pA1/p,  when  the  sign  changes  to  — ,  and  this  seems  to  occur  at  relative  inci¬ 
dences  a/ 0C  between  0.2  (for  slender  bodies)  and  about  0.4  (for  hypersonic  conditions). 

Now  consider  equation  (9d')  at  the  plane  of  symmetry  and  at  the  cone  surface 


Gu  +  £ 


,v]  -  -  [u 


V  +  w„ 


and  substitute  for  w0  from  (C-2) 


ir|>u  +  *rV]  =i[  Wv2~Wp 


v  —  2  u  G 


Now  u  and  v  are  always  positive  on  the  cone  surface  in  the  sign  convention  used,  hence,  when 

d 

the  +  sign  in  (C-2)  applies,  it  follows  that  ~  (£z  u  +  £r  v)  <  '0  and  the  singularity  stays  on  the  surface. 

3  t’5  .  ■ 

When  the  —  sign  in  (C-2)  applies,  —  (£z  u  +  £r  v)  is  still  negative  until  an  incidence  is  reached  such  that 


v2  -  4  p 00/p  >  v  +  2  u  G 


=  3  v 
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since  on  the  surface  at  the  plane  of  symmetry  from  (6)  we  have  v  =  uG.  Alternatively,  this  condition 
can  be  written 

P00  <  -  2  Pv2 

or  in  terms  of  surface  Mach  number  M0,  surface  pressure  p„  and  cone  semi-angle  0<: 

P00  <  -  27  Pe  Ml  sin2 

This  is  the  expression  formed  in  Reference  9  for  the  singularity  to  move  away  from 
£  =  0,  4>  =  t. 

C.1.2  To  Show  That  the  Singularities  Cannot  Move  from  the  Planes  of  Symmetry. 

When  the  singularity  does  move,  and,  in  fact,  for  incidences  above  that  incidence  when 
the  sign  in  (C-2)  first  changes  (i.e.  v2  =  4pfl(j/p  at  the  leeward  generator),  which  always  occurs 
before  adverse  pressure  starts,  we  know  that  the  —  sign  applies  in  the  expression  (C-2)  for  w(<  at  the 
leeward  plane,  that  is 

W0  =  i[-v-V/v2-  4  p00/p  ] 

which  is  always  negative.  On  the  other  hand,  at  the  windward  plane  we  have 

W0  =  i  [  -  v  +  Vv2  ~  4p 00/p  ] 


and  it  is  observed  from  results  that  a  pressure  maximum  always  exists  at  <j>  —  0,  so  that  pp(i  <  0, 
thus  w0  >  0.  Hence,  it  is  shown  that  at  the  surface  w„  >  0  on  the  windward  plane  and  w,,  <  0  on 
the  leeward  plane,  and  since  w  is  zero  at  <f>  =  0  and  ir,  other  zeros  in  w  could  only  occur  on  the  surface 
if  they  were  to  occur  in  pairs.  Now  a  zero  in  w  can  only  occur  at  a  surface  pressure  minimum  or 
maximum,  and  apart  from  the  minimum  and  maximum  pressure  at  the  windward  and  leeward 
generators,  it  is  observed  from  results  that  only  one  other  minimum  occurs  at  incidences  greater 
than  that  at  which  adverse  pressure  first  starts.  Hence  w  cannot  be  zero  except  at  the  windward  and 
leeward  generators. 


Thus  the  singularity  situation  for  a  circular  cone  at  incidence  is  as  follows.  A  saddle  type 
singularity  is  always  situated  at  the  windward  generator  and  a  nodal  singularity  is  at  the  leeward 
generator  until  an  incidence  is  reached  when  p00  <  —  2 7  pi;  M0  sin”  6a.  For  incidences  above  this, 
the  nodal  singularity  will  be  situated  off  the  surface  along  the  leeward  plane  of  symmetry. 


These  conclusions  on  the  position  of  the  singularity  are,  in  fact,  verified  by  results  obtained 
in  certain  cases  for  the  circular  cone  at  incidence,  and  a  plot  of  ^  (£,u  +  £r  v) 

12.5°.  Figure  16  shows  a  plot  of  p,)(4 
a/6c  for  the  same  case. 


and| 


IrV) 


against  £  along  <f>  =  tt  is  given  in  Figure  15  for  the  case  MOT  =  1.797,  d0  = 
2p  v2  at  the  leeward  generator  against  relative  incidence 


It  should  be  noted  that  it  is  possible  by  the  present  method  to  obtain  solutions  for  some 
cases  where  the  singularity  moves  off  the  surface,  since  extrapolation  of  £z  u  +  £r  v  +  from 

£  =  5£  to  £  =  0,  as  described  in  Section  3.0,  is  made.  It  also  appears  to  be  valid  to  extrapolate  press¬ 
ure  through  this  singularity.  The  singularity  does  not  cause  trouble  until  it  moves  so  far  out  that  it 
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is  at  a  distance  corresponding  to  £  approaching  <$£  from  the  surface,  so  that  results  can  still  be 
obtained  for  moderately  large  relative  incidences  (assuming  that  u  +  £r  v  and  p  are  continuous 
through  these  singularities) . 


C.2.0  ELLIPTIC  CONE  AT  ZERO  INCIDENCE 


Since  the  elliptic  cone  at  zero  incidence  is  symmetric  about  <j>  =  0  and  </>  =  ir/2,  only  the 
section  0  <  4>  <  ir/2  of  the  flow  field  needs  to  be  considered. 


Consider  the  pressure  distributions  around  elliptic  cones  in  the  section  0  ^  4,  <  n/2. 
They  have  the  forms  shown,  for  example,  in  Figure  17,  which  shows  the  form  variation  against 
various  values  of  a/b  where  a  is  the  semi-major  axis  and  b  is  the  semi-minor  axis.  The  minor  axis 
lies  in  the  plane  <f>  —  0,  ir.  From  these  distributions  it  can  be  seen  that  the  pressure  is  monotonic 
except  in  the  cases  a/b  =  2.4  and  2.8.  It  is  the  opinion  of  the  author  that  the  non-monotonic  behaviour 
in  these  latter  cases  is  due  to  “numerical  overshoot”  of  the  pressure,  since  the  numerical  procedure  is 
trying  to  fit  an  almost  constant  pressure  in  the  region  0  <  <f>  <  40°  followed  by  a  rapidly  varying 
pressure  in  the  region  40°  <  <f>  <  90°  —  in  fact,  in  the  author’s  opinion,  the  pressure  will  always  be 
monotonic  for  all  a/b  whatever  the  free  stream  Mach  number. 


On  the  assumption  that  the  pressure  is  monotonic  in  the  region  0  <  <t>  <  t/2,  it  is  easily 
seen  that  w  cannot  be  zero  except  at  <p  =  0,  ir/2  as  follows.  If  w  =  0  on  the  surface,  then  from 
equation  (9b)  we  must  have  ^  =  0,  and  hence  from  equation  (9c)  a  necessary  condition  for  w  to 


be  zero  is  that  — ^  must  be  zero.  However,  for  the  elliptic  cone  at  zero  incidence,  since  p  is  monotonic, 

0(p 

w  can  only  be  zero  at  0  —  0,  7r/2. 


The  above  analysis  contradicts  Reference  22,  where  results  are  given  that  show  that  for 
Mw  =  10,  a  =  0.7  and  a/b  =  2,  a  nodal  type  singularity  existed  on  the  cone  surface  at  about 
<t>  =  20°  with  saddle  type  singularities  at.  0  =  0  and  vr/2.  Plots  of  some  results  for  this  same  case, 
using  the  present  method,  are  given  in  Figure  18,  They  indicate  a  nodal  type  singularity  at  <f>  =  0 
with  a  saddle  type  singularity  at  0  =  tt/2. 


It  still  must  be  determined  whether  the  singularities  can  move  off  the  surface  for  the 

elliptic  cone  at  zero  incidence.  It  can  be  shown  that  for  the  elliptic  cone  the  formulae  for  wfi  and  — - 

Of 

(£zu  +  £r  v)  given  by  (C-2)  and  (C-3)  apply  at  the  planes  of  symmetry  at  the  cone  surface,  but  the 
appropriate  signs  have  to  be  determined.  Clearly  when  a/b  =  1,  p00  =  w0  =  0  everywhere  and, 
hence,  the  +  sign  holds  in  (C-2)  and  the  —  sign  in  (C-3).  Now  if  we  assume  a  sign  change  takes 
place  continuously  in  w0  as  a/b  increases, .  this  would  imply  that  v2  —  4p 00/p  =•  0  would  be  a 
necessary  condition  for  the  change,  but  at  0  =  r/2  it  is  observed  that  p,w  <  0  always  for  a/b  >  1, 
hence  v2  —  4p00/p  is  never  zero,  implying  that  the  —  sign  is  always  applicable  in  equation 

(C-3),  This  in  turn  implies  that— ^  u  +  £r  v^J  ^  Q  ^  /2  is  always  negative.  Also,  since 

u  +  £r  vj  ^  is  zero  on  the  cone  surface  and  is  negative  at  the  shock,  it- follows  that  this 
singularity  stays  on  the  surface.  On  the  other  hand,  at  0  =  0,  since  p00  >  0  always  for  a/b  >  1,  then 
clearly  whichever  sign  applies  in  equation  (C-3)  it  follows  that  -(£„  u  +  £r  v) 

<  0,  and  by  the  same  argument  as  above,  the  singularity  stays  on  the  surface. 


{  =  0,  0  =  0 


Thus,  in  conclusion,  the  singularities  stay  at  the  cone  generators  in  the  planes  of  symmetry 
for  the  elliptic  cone  at  zero  incidence. 
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APPENDIX  D 

THE  SOLUTION  FOR  A  CIRCULAR  CONE  AT  ZERO  INCIDENCE 


The  computations  to  generate  the  zero  order  solution  are  straightforward  once  a  reasonable 
initial  estimate  of  the  shock  angle  is  determined.  The  method  will  be  described  without  going  into 
mathematical  details. 


An  estimate  of  the  shock  angle  is  made  empirically  from  results  obtained  by  Sims  (Ref.  2) 
for  situations  involving  conically  slender  bodies  (M^  sin  0„  <  0.6)  and  from  hypersonic  small 
disturbance  theory  for  Mro  sin  0Q  >  0.6.  These  estimates  then  take  the  form 


Mod  sin  0a 


-  1  = 


y+1 

2 


(Mod  sin  0„)3'63 


for  Mco  sin  0C  <  0.6 


and 


sin  0  r7  +  1  1 

sin  0C  2  M 3  sin2  0C 

CO 


where  0S  is  the  shock  wave  angle. 


for  Mco  sin  0C  >  0.6 


Having  obtained  these  estimates  of  the  shock  angle,  quantities  behind  the  shock  can  be 
determined  and  the  equations  of  motion  integrated  into  the  body,  thus  giving  a  residual  normal 
velocity  on  the  surface.  An  iteration  procedure  due  to  Wegstein,  (Ref.  23)  is  then  used  to  improve 
the  estimate  of  the  shock  angle  so  that  the  residual  normal  velocity  is  made  very  close  to  its  required 
value  of  zero.  It  is  found  that  this  computation  requires  only  three  or  four  iterations,  hence  only  a 
few  seconds  of  computer  time. 
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APPENDIX  E 

SOME  COMMENTS  ON  THE  COMPUTER  PROGRAM 


E.1.0  THE  INCLUSION  OF  FURTHER  FOURIER  COEFFICIENTS  TO  REPRESENT  THE  SHOCK  SHAPE 


As  the  incidence  or  other  perturbation  is  increased,  the  number  of  Fourier  coefficients  in 

m 

the  shock  wave  equation  F  (0)  =  22  Fi  cos  i  0  must  be  increased.  At  each  incidence  the  estimated 

•  =  0  m 

Fm  (i.e.  the  last  coefficient)  is  compared  with  22  Fi  (the  tangent  of  the  shock  angle  at  0  =  0), 


i  =  0 


and  if  |  Fm  1/22  F(  >  5  xTO'6,  an  extra  coefficient  is  added  so  that  the  shock  wave  equation  is 


m  +  1 

now  F  (0)  =  22  cos  i  4>  where  F0,  Fi . Fm  have  the  same  values  as  were  estimated  previously 

i  =  0 

and  the  estimate  for  Fm+  1  is  zero.  An  extra  coffieeient  is  also  added  if  the  sum  of  squares  of  residuals 
on  the  preceding  calculation  is  greater  than  2  x  10'6. 


E.2.0  STEP  LENGTHS 

Step  lengths  of  22|°  and  0.1  in  </> and  £  respectively  seem  to  give  sufficiently  accurate  results 
for  all  the  cases  computed  so  far  for  the  circular  cone  at  incidence.  To  keep  the  step  length  of  22|° 
in  4>  implies  that  no  more  than  nine  Fourier  coefficients  can  be  used  to  represent  the  shock  shape, 
since  the  number  of  unknowns  defining  the  shock  must  not  exceed  the  number  of  dividing  lines  in 
the  region  0  <  0.<  w/2.  Thus,  when  the  criterion  in  Section  E-1.0  is  satisfied  to  include  a  tenth 
coefficient,  the  Fourier  series  representation  to  the  shock  is  not  used;  instead  the  shock  is  simply 

represented  by  values  at  discrete  points  0=0,  22i . 180°.  Initial  estimates  at  these  discrete 

points  are  found  by  extrapolation  as  before.  The  iteration  procedure  then  improves  these  estimates 
at  the  discrete  points. 

In  the  cases  of  the  elliptic  cones  and  the  fourth  order  cosine  Fourier  series  body,  the 
increment  in  4>  is  taken  as  22.5°  to  start  with,  i.e.,  when  the  configuration  is  near  to  a  circular  cone 
at  zero  incidence  (this  is  equivalent  to  8  sections  and  9  lines  in  the  range  0  <  0  <  180°).  As  the 
incidence  or  eccentricity  or  other  perturbation  increases,  further  Fourier  coefficients  are  added  to 
the  sum  forming  the  shock  r  =  F(0),  and  when  the  number  of  coefficients  is  about  to  exceed  the  num¬ 
ber  of  dividing  lines  in  the  region  0  <  0  ^  180,  the  number  of  lines  is  increased  so  that  their  number 
is  always  equal  to  the  number  of  coefficients.  It  is  necessary  to  do  this  in  order  that  the  number  of 
discrete  normal  velocities  at  the  surface  is  not  less  than  number  of  unknowns  defining  the  shock, 
which  is  a  necessary  condition  for  the  iteration  process. 
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E.3.0  EXTRAPOLATION  TO  THE  BODY 

Integration  was  always  made  to  a  value  £  =  0.1  in  equal  steps  5f.  The  velocity  u  f  fr  v 
+  -^w  was  then  extrapolated,  using  the  values  and  derivatives  of  this  quantity  at  0.1,  0.1  +  of, 
0.1  4-  2  5f,  and  0.1  +  3  <$£,  by  solving  the  formulae 

q(0)  =  q  (Sh)  -  5h  q'  (ah)  +  [q"(0)  +  ah  q"'(0)  +  q-(0)  J 

ah3  r  n  tu-t 

-  [  q"'(0)  +  ah  q*v(0)  J  +  -g*  q-v(0)  +  0  (5h6) 

where  6h  =  0.1  —  (j  —  1)  5f  for  j  =  1,  2,  3,  4  (af  <  0) 

and  q  represents  the  quantity  to  be  extrapolated. 

When  the  iteration  procedure  is  completed  the  differential  equations  are  integrated  to  a 
value  f  of  approximately  0.003  and  the  pressure  is  then  extrapolated  Vo  the  surface  by  a  quadratic 
formula. 

E.4.0  COMPUTING  TIMES 

Typical  times  on  an  IBM  360/50  computer  for  the  circular  cone  at  one  incidence  are  40 
seconds  to  1  minute  for  relative  incidences  up  to  0.6;  times  for  higher  relative  incidences  vary  from 
1  to  3  minutes  approximately.  Actual  times  are  given  below  for  certain  typical  cases. 


\ 
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E.4.1  Circular  Cone 


@0  =  25° 


e a 
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a 

£■ 
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E.4.2  Elliptic  Cone 


Mco  =  10 


a  =  0.7 


a/b 


Time  (min.) 


1.01 

1.2 

1.4 
1.6 
1.8 
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2.8 


0.5 
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1.2 
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b  =  0.18475 
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E.5.0  TERMINATION  OF  THE  ITERATION  PROCEDURE 

The  iteration  procedure  is  continued  until  the  change  in  each  Fourier  coefficient  is  less 
than  10 "b  of  that  Fourier  coefficient.  An  alternative  criterion  for  terminating  the  iteration  is  also 
used;  this  is  based  on  the  sum  of  squares  of  residual  errors  being  less  than  10  b.  It  was  found  that  for 
small  relative  incidences  up  to  about  0.6,  only  one  step  in  the  iteration  procedure  was  required  to 
give  sufficient  accuracy.  This  one  step  requires  about  10  to  15  integrations  of  the  equations  from  the 
shock  to  the  body.  For  higher  relative  incidences  about  2  to  4  steps  were  necessary,  requiring  about 
15  to  30  integrations  of  the  equations. 

E.6.0  PROGRAM  AVAILABILITY 

Duplicate  program  decks  are  available  upon  application  to  the  High  Speed  Aerodynamics 
Section,  National  Aeronautical  Establishment. 

A  description  of  the  program  and  its  use  is  given  in  Reference  24. 


